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ABSTRACT 



We have developed a relativistic, radiation-hydrodynamics Lagrangian code, specificaUy tailored to 
simulate the evolution of the main observables (light curve, evolution of photospheric velocity and 
temperature) in core-collapse supernova (CC-SN) events. The distinctive features of the code are 
an accurate treatment of radiative transfer coupled to relativistic hydrodynamics, a self-consistent 
treatment of the evolution of the innermost ejecta taking into account the gravitational effects of the 
central compact remnant, and a fully implicit Lagrangian approach to the solution of the coupled 
non-linear finite difference system of equations. Our aim is to use it as numerical tool to perform 
calculations of grid of models to be compared with observation of CC-SNe. In this paper we present 
some testcase simulations and a comparison with observations of SN 1987A, as well as with the 
results obtained with other numerical codes. We also briefly discuss the influence of the main physical 
parameters (ejected mass, progenitor radius, explosion energy, amount of ^^Ni) on the evolution of 
the ejecta, and the implications of our results in connection with the possibility to "standardize" 
hydrogen-rich CC-SNe for using them as candles to measure cosmological distances. 
Subject headings: supernovae: general — hydrodynamics — radiative transfer — methods: numerical 
— supernovae: individual (SN 1987A) — distance scale 



1. INTRODUCTION 

It is widely accepted that core-collapse supernovae 
(CC-SNe) represent the final explosive evolutionary 
phase of stars having initial (i. e. main sequence) masses 
larger than - 8-10Mf7, (e g. l\yooslev fc Weavej 119861: 
Eamuy 200331 IHeeer et al.ll2003D . As such. CC-SNe are 
fundamental probes of the stellar evolution of massive 
stars and can be used to understand the link among 
explosion mechanism, nature of remnants, progenitors 
evolution, and environment around pro genitors (e.g. 
lFilippenkolfl997l : iCappellaro fc Turatto( l2b00') . 

In addition to their intrinsic interest, CC-SNe are 
relevant to many astrophysical issues. For exam- 
ple, the CC-SNe influence the physical evolution of 
galaxies, their frequency bein g related to the on- 



going star formation rate (e.g . ICappellaro et al.l 120051 : 
iMadau. Delia Valle. fc Panagial I1998D and determining 
the kinemati cs and the structure of the inter stellar 
medium (e.g. IRatnatunga fc vaii den Bergh|[T98l . CC- 
SNe are also important because of their role in the 
production of neutrinos, cosmic rays and, probably , 
gravitation al waves ( e.g. iHaungs. Rebel, fc Rotlil 120031 : 
iDuan fc K ncUer 20p iPagharoUet al.ll2009D . Addition- 
ally, the CC-SN ejecta, enriched in heavy elements, 
make CC-SNe key actors in the nucleosynthesis pro- 
cesses of intermediate and trans-iron ele ments and in 
the chemical evolution of g a laxies (e.g. lArnettI 119961 : 
iChieffi. Limongi. fc Stranierol I1998D . Moreover, they 
seem to be particularly promising to measure cos- 
mological distance s in a ddition to ty pe la SNe (e.g . 
INugent et al.l 120061: [Zainpierill2007t loiivares et al.l 12010 . 
and references therein). 
In spite of the importance of these explosive events 



in astrophysics, there are still basic questions to be an- 
swered on CC-SNe, related to the extreme variety of their 
displays. Indeed they appear to show different energetics 
and to eject rather diverse amounts of material, causing 
rather heter ogeneous behaviors of their light curv es and 
spectra (e.g. iTuratto. Benetti. fc Pastorelloll2007f) . This 
apparent heterogeneity may be linked to stel lar evolution 
effect s and to the explosion mechanism (e.g. lPumo et al.l 
120091 and references therein), but the exact link between 
the physical properties of the explosion (ejected mass, 
explosion energy, stellar structure and composition at 
the explosion) and the observational characteristics is far 
from being well-established. 

Light curve and spectral modelling have often 
successfully been used to provide information about 
the physical properties of single CC-SNe (e.g. SNe 
1987A, 1993J, 1999em, 2004et and 2005cs; see 
Shigevama. Nomoto. fc Hashimotolll988t lWooslevl | 1988 : 

. It 



Shigevama fc Nomotd '1990': 'Blinniko v et al.l 119981: 
Utroibini 12 004. 2007; Baklanov, Blinni kov. fc PaylvukI 



20051: lUtro bin fc Chugai 2008; Pas torello et al.l 120091: 
Bersten. Ben venuto. fc Hamuy 20U)- Comparatively 
little effort has been devoted to investigating large 
samples of these events and trying to explain the signifi- 
cant range of properties that they show. Differences in 
luminosity, expansion velocity and ^®Ni yields are very 
large (up to ^ 2 order of magnitude), and difficult to 
relate to changes in a single parameter, such as the mass 
of the progenitor star. However, despite this, variations 
in the observed properties appear to obey a certain 
order. Some correl ations among differe nt observables 
have been noted (e.g. lHamuvll2003bl ; iZamp ieri 2007) and 
used to calibrate CC-SNe as distance indicators (e.g. 
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Hamuv fc Pintoll2002l: lHainuvll2003bl: iNugent et aL|[200l 



Olivares et alj I2010D . in addition to Baade-Wesselink 
type, spectral based techniq ues (like t he expanding 
photosphere method, see e.g. iKirshner & : Kwan 197; 
Eastman. Schmidt, fc Kirshneri Il996t IDessart et al 



2. CODE DESCRIPTION 

The code is a new, improved version of a general- 
relativistic, radia tion hydrodynamics, Lag rangian 



20081 or the spectral-fitting expa,nding atmospher e 
method, see e.g. iMitchell et al][200a iBaron et al.|[200l . 

While parameter studies of the light curves of CC- 
SNe have been largely pursued in the p a st both with 
analytical (e.g. IBlinnikov fc Popovl 119931: IPopovl 119931: 
[Arnet t 1996, and refer ences therein) a nd numerical 
(e.g. iLitvinova fc N adczhin 1983, 1985'; lYound [200l 
IDessart. Livne. fc \Valdman„201Q,) techniques, only a few 
systematic theoretical investigations of the physical ori- 
gin of the aforementioned correlations have been per- 
formed (e.g. [^ampicri 2007; Kascn & Wooslcv 2009). It 
is still not completely clear how variations of basic pa- 
rameters conspire in such a w ay to give rise to the ob - 
served relations. In particular, iKasen fc WooslevI ()2009f ) 
found that the relation between luminosity and expan- 
sion velocity at 50 days is explained by the behavior 
of hydrogen recombination in the SN envelope and pro- 
posed additional correlations based only on photometric 
data. A better physical understanding of these correla- 
tions would allow us, on one side, to pinpoint possible 
biases that may affect the use of CC-SNe as standard 
candles and, on the other side, to look for other corre- 
lations that may reveal more promising when applied to 
high redshift SNe. 

With the aim of clarifying the aforementioned 
astrophysical issues, we have improved on previ- 
ous work that we started on the numerical mod- 
elling; of SN eiecta and fallback (iZampieri et al.l 119981 : 
iBalberg. Zampieri. fc Sh apiro 200JJ), and on the physi- 
cal interpretation of the h eteroge neous b ehavior of Type 
II pl ateau SNe (jZampierf et all 120031: IZampieril 120051 
I2007D . In this paper we present a new, improved 
version of a general-relativistic, radiation hydrodynam- 
ics, Lagrangian code tailored to the modelling of CC- 
SNe. Similar numerical treatment of the expanding 
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) . Our code has the advantage of being able to fol- 



low the fallback of material on the central remnant in 
a fully relativistic formalism, and allows for an accurate 
treatment of radiative transfer in all regimes. It will en- 
able us to explore the physical behavior of CC-SNe, to 
link their physical and observational properties, to inves- 
tigate the existence of correlations among their different 
observables, and also to perform model fitting of single 
events. 

The plan of the paper is the following. In sect.[2]we de- 
scribe the features of the new version of the code (details 
on the adopted numerical procedure and on the finite 
difference form of the equations are presented in the Ap- 
pendix). Sect. [3] illustrates the numerical calculations. 
Sect, mis devoted to the code validation and the analysis 
of the basic physical properties of the ejecta in presence 
of a central compact remnant during all the different evo- 
lutionary stages. A summary with final comments and a 
short discussion of the implications of our results in con- 
nection with the possibility to "standardizate" hydrogen- 
rich CC-SNe is presented in Sect. [5l 



code described in IZampieri Miller, fc TuroUal (Il996.) . 
IZampieri et al.l ()1998D and in IBalberg et aTT " poooir and 

originally developed for studying spherical accretion 
onto black holes and fall back in the aftermath of a SN 
explosion. 

This new version of the code is able to compute the 
evolution of the ejecta and the emitted luminosity, from 
the breakout of the shock wave at the stellar surface up to 
the nebular stage (when the envelope has recombined and 
the energy budget is dominated by the radioactive decays 
of the heavy elements synthesized in the explosion). We 
solve the equations of relativistic radiation hydrodynam- 
ics in spherically symmetry for a self-gravitating matter 
fluid which interacts with radiation, taking into account 
the heating effects due to decays of radioactive isotopes 
synthesized in the SN explosion. The distinctive features 
of the code are summarized here, while a detailed de- 
scription of the equations, numerical method and input 
physics is reported below (see sections \n\ through I2.4p . 
The code encompasses: 

- an accurate treatment of radiative transfer in all 
regimes from the one in which the ejecta are opti- 
cally thick up to when they are completely trans- 
parent. This is obtained through the solution 
of the first two moments of the radiative trans- 
fer equations (Eqs. [5] and [3] in sect. 12.11) . Al- 
though it is necessary to adopt a closure rela- 
tion, which is constrained only in the asymptotic 
optically th ick /thin regimes (see IZampieri et al.l 
IT9961 IT998L for more details), the uncertainty in- 
duced by this approximation is not significant for 
wavelength( frequency)-integrated ra diative trans- 
fer (see e.g. iNobih fc Turollal I1988D and this ap- 
proach is definitely superior to any treatment based 
on the diffusion approximation; 

- the coupling of the radiation moment equations 
with the equations of relativistic hydrodynamics, 
adopting a fully implicit Lagrangian finite differ- 
ence scheme (see sect. 12. ll and 12.2! for details); 

- the possibility to compute the evolution of the 
ejecta and the emitted luminosity taking into ac- 
count the gravitational effects of the compact rem- 
nant. This enables the code to deal with the 
fallback of material onto the compact remnant 
and, consequently, to accurately determining the 
amount of ejected ^^Ni. 

The main limit ation of our code compared to other 
simi lar codes (e.g. lYound 12004 IKasen fc WooslevI I2009L 



and IDessart, Livne. fc WaldmanI I2010D is that we con- 
sider only wavelength integrated radiation energy density 
and luminosity, and hence we cannot compute the ac- 
tual spectrum of the SN. Also, the evolution starts from 
"ad hoc" initial conditions that mimic the profiles of the 
physical variables after shock and reverse shock passage, 
but such conditions are not meant to be "perfectly realis- 
tic" , in the sense to be the outcome of a stellar evolution 
code coupled with numerical computations of the explo- 
sion phase. Some improvements in this respect have been 
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already added to the present version of the code (see 
sect. 12.31 for details). At the same time, we are work- 
ing (in collaboration with A. ChiefR and M. Limongi) 
at interfacing it with the output of a hydrodynamical 
code that follows the CC-SN explosion and the explosive 
nucleosynthesis, starting from pre-SN models obtained 
through a stellar evolutionary code. 

2.1. Equations 

Since we consider the effects of the central compact 
remnant on the evolution of the ejecta, we adopt a fully 
general relativistic treatment. The equations govern- 
ing the relativistic radiation hydrodynamics of the ex- 
panding ejecta i n sphericall y sym metry are reported i n 
iZampieri etUI (|1996L 119981 ) and iBalberg et all (|2000l ). 
The main physical variables that describe the behavior 
of the ejecta are the gas mass density p, and the gas in- 
ternal energy per unit mass e and pressure p (measured 
in the frame comoving with the ejecta). Those describ- 
ing the radiation field are the radiation energy density 
W o and flux wi (both i n units of erg cm~^). We refer 
to IZampieri et al.l (|1998l ) for the full set of equations and 
related quantities, while here we summarize for conve- 
nience the gas energy equation and the first two moments 
of the radiative transfer equations, whose treatment has 
been deeply modified in the present analysis. They read: 

e,t+akp{B ~wo)+p[-] = Q energy eq., (1) 
\Pj,t 



(wo),t - akpp{B - Wo) + 
4/6* t\ , fbt rt 

3 T+^y +nT-T 



Wo 

1 



(wia^r^).p 



zero-th moment eq., 



{wi)^t + akftpwi + 



(2) 



2z.r(> + ^ 
r 



I 



■{woa'^),f,+ 



1 



35r3 



{fwoar^).f. 



3a% 

first moment eq. 



(3) 



where t and p are the Lagrangian time and mass con- 
tained within a comoving spherical shell of radius r (a 
comma denotes partial derivates with respect to the cor- 
responding variable), a (00 metric coefficien t) is a func- 
tion of t and p (computed from Eq. [6] in Za mpieri et al.l 
[T998I). 6 = l/(47rr2p), B = aflT* (an blackbody radia- 
tion constant), Q is the heating rate of radioactive decays 
from the isotopes synthesized in the SN explosion (see 
sect. 12. 4p . and kp and kp are Planck mean and Rosse- 
land mean opacity, respectively. The dependence on the 
gas temperature T is contained in e, p, kp and kp, and 
is specified through the equations of state (see sect. 12. 4p . 
kp and kp are calculated interpolating the TOPS opaci- 
ties as a function of T and p (see again sect. 12. 4p . Finally, 
the function / (Eddington factor), relating the second- 
order moment of the radiation intensity to wp, i s calc u- 
lated using Eqs. [12] and [13] in , Zampieri et al.l (|1998l ). 



2.2. Numerical method 

In order to solve the complete system of equations of 
relativistic radiation hydrodynamics, the old version of 
the code adopts a semi-implicit Lagrangian finite differ- 
ence scheme where the time step is controlled by the 
Courant condition and the requirement that the frac- 
tional variation of the variables in one time-step be 
smaller than 10%. The energy equation (Eq. [Tj) and 
the zero-th moment of the radiative transfer equation 
(Eq. [5]) form a non-linear system of equations in the 
gas temperature T and are then solved point-by-point 
on the computational grid using a Newton-Raphson iter- 
ative method (see the Appendix in IZampieri et al1ll998L 
for details). 

We have deeply modified the numerical treatment pre- 
viously adopted by implementing a fully implicit La- 
grangian finite difference scheme. This allows for a major 
improvement in the numerical stability and overall com- 
putational efficiency of the code, especially during those 
phases when fast motions of steep gradients occur (e.g. 
the radiative recombination phase). The first moment 
equation for the radiative flux wi (Eq. [3] ) is now solved 
together with Eqs. [Tj and [5] in a fully implicit scheme on 
the whole computational grid, that requires a modifica- 
tion of the temporal centering of the variables. This leads 
to a highly non-linear system of equations that is solved 
through a Newton-Raphson iterative method and matrix 
inversion packages that minimise the CPU time and the 
required storage space. In the Appendix we report more 
details on the finite difference form of the Eqs. [1], [2] 
and [3], and the numerical procedure adopted to solve 
them. 

2.3. Initial and boundary conditions 

In the present analysis we consider the evolution of the 
ejecta after the explosion phascQ. At this stage, the SN 
shock wave has already propagated through the envelope 
of the progenitor star redistributing the explosion energy 
through it, and the evolution starts when the envelope 
is essentially free-coasting and in homologous expansion. 
We assume that it is comprised of a mixture of hydrogen, 
helium and heavier elements. As oxygen is expected to 
be the most abundant heavy element, we further assume 
that it represents the entire m etal component in t he en- 
velope's final composition (see IBalberg et ani2000l ). The 
mass fraction of ^^Ni is assumed to be concentrated to- 
wards the center and to vary as a function of Lagrangian 
mass p as: 



X56Ni(M) = ^^<>Ni,o X exp 



(4) 



where Xsep^j q is the initial mass fraction of ^®Ni at the in- 
ner boundary of the ejecta, Kmix is a numerical constant 
whose value is essentially related to the extent of mixing 

^ We are working at interfacing our code with the output of 
self-consistent calculations of the explosion phase, but during this 
testing phase we preferred to assume idealized initial conditions 
that provide an approximate description of the ejected material 
after shock (and possible reverse shock) passage. Clearly, as the 
actual velocity, density and heavy element distributions of the post- 
shock material may affect the light curve and the determination of 
the envelope's physical parameters in a significant way, at present 
we can provide only approximate estimates of these quantities. 
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of ^^Ni throughout the envelope, and ^imax is the value 
of /i at the outer boundary, equal to the total envelope 
mass (see Appendix |X] for details) . 

Assuming that, at the onset of expansion, radiation 
is in local thermodynamic equilibrium (LTE hereafter) 
with the gas throughout the envelope and the initial pho- 
ton mean free path Xq is much smaller than the radius of 
the envelope, the initial temperature profile in the ejecta 
can be well approxiniated by the so-called "radiative zero 
solution" of lArnettl (|1980( ): 



sm(7rx) 



(5) 



where Tq is the initial central temperature, x = r/Ro 
and i?o = ^{p-max) is the initial radius of the outermost 
shell of the ejecta. 

As the gas is initially in LTE with radiation, we assume 
that a polytropic relation with index F = 4/3 holds. For 
a radiation dominated flow, P = anT^/i — Kp^l^ {K 
polytropic constant) and hence 



p{r{fi),t = 0) = po X 



sin{'Kx) 



nx 



(6) 



with pQ = (afl,/3i^)3/4r3. 

Considering that the envelope is expanding homolo- 
gously, the initial velocity profile is given by: 



j{r{p),t = 0) = Vo X 



Rr, 



(7) 



where Vq is the initial velocity at the outer radial bound- 
ary Rq. 

As for the radiation field, in the aforementioned as- 
sumptions the diffusion approximation holds. There- 
fore, wq — urT'^ and wi « wq/ (knpr), where kn is ap- 
proximated with the electron sca ttering opacity k^s (see 
iNobih. Turolla fc ZampiedHlOQl . 

For a fixed composition and radial distribution of 
^^Ni, four input parameters determine uniquely the dy- 
namical and thermal properties of the expanding CC- 
SN envelope (including its initial kinetic and ther- 
mal energy): the total mass Me„„, the initial ra- 
dius _Ro: the initial sound speed at the inner bound- 
, and the ratio k of the initial accretion 
timescale tacc,o to the initial expansion time tea^p.n , 
defined as foll ows (IColpi. Shapiro. &: WassermanI 119961 : 
iZampieri et"allll998l )l 



'^acc,0 



GAL 



Rq 



(8) 



where Ale is the mass of the compact remnant, taken to 
be equal to 3Mq at the onset of the evolutiorQ. In par- 
ticular, the polytropic constant K and the initial central 

^ This value of is larger than that of a typical com- 
pact remnantensuing_fr CC-SN events {~ 1.3-2.6M0; see 
e.g. ILimongi. Straniero fc Ch ieffi 2000), but we chose it because 
it helps the numerical stability of the innermost part of the flow. 
While the dynamical evolution of the expanding ejecta is essen- 
tially unaffected by Mc, the amount of material that falls back 
in the compact remnant during the evolution is slightly overesti- 
mated by a factor < 3-8%. Clearly, detailed calculations of fallback 
require choosing the value of more accurately. 



temperature Tq depend on the three input parameters 
Cs„, Rq and Mem,, and can be written as K — 30^^/4^^^ 

and To — (SK/ajij^^^pQ^^ , where pq is a function of Rq 
and A'/enti- Vq can be expressed in terms of Cgg, Rq and 
k through Eq. [5]. The dynamical interplay of the in- 
ner accreting envelope and the outer expanding ejecta 
can be characterized, in absence of ^^Ni energy input, 
in terms of the sole parameter k (see IColpi et all 119961 : 
IZampieri et al.|[l998l for details), while the initial kinetic 
and thermal energy of the envelope are obtained by in- 
tegration. 

Concerning the boundary cond itions and their numer- 
ical implementation, we refer to IZampieri et all ([l998). 
However we recall that at the inner boundary we assume 
negligible pressure and radiative forces on the gas, and 
adopt an outgoing wave (floating) boundary for wi. At 
every time step we also set the mass of the compact rem- 
nant equal to its initial mass plus the actual mass that 
has fallen back up to that time. At the outer boundary, 
we adopt free expansion for the ejecta, a non-illuminated 
atmosphere in radiative equilibrium for the radiation 
field, and synchronization of the coordinate time with 
the proper time of an observer comoving with the ejecta. 

2.4. Input physics 

While a complete description of the input physic s 
adopted in the paper is reported in lBalberg et al.l (|200C1( ). 
here we summarize the most relevant aspects concerning 
the treatment of the Rosseland and Planck opacities, the 
ionization fractions, the gas equations of state, and the 
radioactive heating. 

As for the opacity and ionization fractions, the code 
uses the TOPS opacities availab le at the LANL s erver 
(version preceding April 2000; see lMagee et al.iri995| ). ex- 
tended in the low- temperature regime (T < 5.8x 10^ K) 
using the tables of I Alexander fc Fergusonl (|1994l ). They 
are calculated with linear interpolation in temperature, 
density and chemical composition, considering oxygen 
as representative of the entire metal component in the 
ejecta. A more refined approach using the most up- 
to-date opacities and realistic post-explosive chemical 
composition (computed following the progenitor evolu- 
tion from the main-sequence up to the explosion) will 
be presented in a forthcoming paper. It should be noted 
that these opacities are computed assuming LTE between 
matter and radiation, while in the SN ejecta these condi- 
tions are not strictly met near to the photosphere. Fur- 
thermore, the effect of the non-thermal excitation and 
ionization from gamma rays is not included in the cal- 
culation of the opacities. In order to account for it, a 
so-called "opacity floor" (i.e. a minimum value equal 
to 0.1 cm^ g~^) is used in all the simulations reported 
in this paper, following a common practi ce in the hy- 
drodynamical mod elling of CC-SNe (e.g. iHe rzig et al.l 
1990; Swartz, Wh eeler, fc HarknessI Il99ll : l^ung 200i 
iDessart et al.ll2010t ). 

Concerning the gas equation of state, it is approxi- 
mated as an ideal gas of ions and electrons. 

As far as radioactive heating is concerned, we assume 
that the energy released by the decays of all the radioac- 
tive isotopes present in a given mass shell is absorbed 
locally or escapes from the envelope. Non-local effects in- 
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TABLE 1 

Properties of the relevant radioactive isotopes 



Isotope 



[erg i 



[erg g s 



[days] 



^^Ni 3.90x10"' 8.8 

s^Co 6.40x10^ 2.24x10* 111.3 

S'^Co 6.81x10* 391.0 

**Ti 2.06x10* 6.54x10^^ 3.28x10"' 

N ote. — Data taken fromlWooslev. Pinto, fc HartmannI II1989I ') 
and lShigevama fc Nomotd 1 11990) . 



duced by gamma-ray heating are not considered. This is 
implemented adding to the gas energy equation (Eq. [Tj ) 
a local energy source term Q, which represents the ra- 
dioactive heating rate of all the relevant isotopes synthe- 
sized in the CC-SN explosion and is given by: 



(9) 



where Xi, Ti, Si^-y and e+ represent the mass fraction, 
lifetime, energy per unit mass and unit time released by 
the decay of the i-th isotope in the form of 7-rays and 
positrons, respectively, and the factor fi accounts for the 
fact that the 7-rays are not totally trapp ed in the enve- 
lope (for details, see lBalberg et al.l 12000 . and references 
therein). The relevant radioactive nuclei considered in 
the present investigation and the values of their char- 
acteristic parameters are reported in Table [T] The ^^Co 
abundance is evaluated considering that the isotope ^^Co 
is involved into the nuclear decay chain ^^Ni — > '^^Co — 
^^Fe, so its mass as a function of the time t is determined 
from: 



Ni.O ^ 



M56Coit) = A/56Co,0 X e ^=«Co + 

1 



1 



e Co 



MseNi.o X e "^•'co - e 



(10) 



where we assume that the initial abundance of ^^Co is 
equal to zero and that 1/(1 — TseNi/rseco) ~ 1. 

3. NUMERICAL COMPUTATIONS 

3.1. Preliminary tests 

In order to test the dependability of the new version 
of the code, a few simulations have been performed with 
both the new and old version. The evolution starts from 
the same set of initial conditions reported in sect. 12.31 ex- 
cept for the radial distribution of ^^Ni that is assumed to 
be uniform. We varied several input parameters (num- 
ber of points, location of the inner boundary, mass of 
the ejecta, initial radius and energy) to check the sta- 
bility and accuracy of the code. In the outer expand- 
ing part of the flow (which comprises most of the mass, 
in practice all the envelope mass apart from the inner- 
most ^ 10~^Mq) the fractional difference between the 
values of the variables computed with the new and old 
version of the code is < 10-20% (e.g.. Fi gures [T] and [H 
see also lPumo. Zampieri fc Turattoll20iol ). The only ex- 
ception is the profile of the radiative flux that differs by 
< 50%. Indeed, the radiative flux depends sensitively on 



the numerical treatment. The improved stability of the 
new version makes the profile computed by the new code 
more reliable. A significant difference (< 40%) is present 
also in the innermost ^ 10"'^ Mq of the envelope (below 
logi? ^ 12.7), where the flow starts to fall back onto 
the remnant. This is most likely linked to the delicate 
balance between the pressure gradients and the gravita- 
tional force, th at determines the lo cation of the accretion 
radius ra (see iBalberg et al.l I2000L for details) . Even a 
slight difference in the calculation of the gas and radia- 
tion pressures may cause a sign reversal in the velocity of 
the marginally bound gas shells (containing < 10~^Mq ) 
located near to r^, that start then to fall back onto the 
remnant instead of going outwards. 




12.5 13 13.5 

log R [cm] 



14 14.5 



Fig. 1. — Comparison between results obtained from the old (dot- 
ted line) and new (dashed line) versions of the code. The compar- 
ison refers to a model evolved for 15 days from the breakout of 
the shock wave at the stellar surface, and having initial radius of 
3 X lO'^^ cm, total energy of 1 foe, amount of ^°Ni equal to 0.07 
Mq, and the envelope mass of 16 Mq. Top, middle, and bot- 
tom panels show the radial profiles of the gas matter density p 
(in units of g cm~'^), absolute value of the fluid radial-velocity (in 
units of the velocity of light c) and the gas temperature (in units 
of Kelvin), respectively. The innermost part of the envelope, below 
log/? ^ 12.5-12.7, is accreting onto the central remnant and has 
negative velocity. 





Fig. 2. — Same as Fig.[T] but for the radial profiles of the radiation 
energy density wo (top panel) and the radiative flux wi (bottom 
panel). Both quantities are in units of erg cm~'^. 
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Fig. 3. — Radial profiles of the radiative luminosity at the begin- 
ning of the recombination phase for model (1) in TableO Less than 
5% in mass of the envelope has recombined. The sharp front at 
logi? ~ 15.1 marks the position of the recombination wave, where 
also the photosphere is located. 

While the simulations with the old version often de- 
velop numerical instabilities, those evolved with the new 
version show to be stable. Several models were evolved 
for long times checking the long-term numerical stabil- 
ity of the new code during some key evolutionary stages 
like the radiative recombination phase. The position 
of the recombination front is very well traced (see Fig- 
ure[31) and, consequently, it is possible to follow the post- 
explosion evolution of a CC-SN event in its "entirety", 
from the breakout of the shock wave at the stellar surface 
up to the nebular stage. 

3.2. Grid of models 

Different simulations were performed with the new ver- 
sion of the code exploring a relatively large fraction of the 
parameter space related to the post-explosion evolution 
of a "typical" CC-SN. 

In particular a total of 22 models were calculated 
changing the initial radius from 3 x 10^^ to 10^*^ cm, 
the total energy from 0.5 to 3 foe, the amount of ^^Ni 
from 0.005 to 0.070 M©, and the envelope mass from 8 
to 18 M0(see Table[2]for details). In the present sample, 
we assume that the initial thermal and kinetic energy are 
equal and hence three out of the four input parameters 
of the model {M^nv, Csq and k) are related. All the sim- 
ulations were performed fixing the parameter controlling 
the ^^Ni distribution at the value Kmix = 5. This cor- 
responds to having ^ 90% of the ^^Ni in the envelope 
concentrated in the innermost ~ 3, 5, 6.5 and 7.5AfQ for 
the models with an envelope mass equal to 8, 12, 16 and 
I8M0, respectively (see Eq. [4]). 

4. VALIDATION OF THE NEW CODE AND 
PHYSICAL BEHAVIOR OF CC-SN MODELS 

4.1. Code validation 

Some of the models reported in Table [2] have been used 
to further validate the code against observations and by 
comparison with similar models reported in the litera- 
ture. In particular, we have compared: 

- model (1) with model (A) of lYound (|2004D that 
has a similar input parameters (i.e. total energy 
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of 1 foe, initial radius of 3 x 10"'^^ cm, ^^Ni mass 
of O.O7M0, and ^^Ni mixing throughout the inner- 
most 6 Mq region of the envelope); 

- models (3), (4) and (10) with results of simu- 
lations pe rformed with the semi-analytic code of 
iZampieri et al. (20Q3) andfZampieri (2007); 

- model (13) with results of simulations performed 
with the aforementioned semi- analytic code, and 
with model (B) of lYounj (|2004l ) that has a similar 
input physics (i.e. total energy of 1 foe, initial ra- 
dius of 3 X 10^3 56Ni mass of O.O7M0, and ^^Ni 
mixing throughout the innermost 6 Mq region of 
the envelope); 

- models (1) and (2) with observations of SN 1987A. 

Figures l4l through [TT] show the results of these compar- 
isons. 

The agreement between the results obtained with our 
code and those reported in the literature is good, consid- 
ering the different initial conditions, input physics and 
numerical methods adopted (see Figures H] through ^ . 
In particular, we found overall c onsistency with t he re - 
sult s of the semi-an alytic code of IZampieri et al.l ()2003f ) 
and IZam pierll (j2007D , that serves also as cross validation 
of the latter. 

However, in some simulations there are somewhat sig- 
nificant discrepancies in the photospheric velocity and 
temperature during the early evolution {< 10-20 days). 
The fractional difference between the values of the vari- 
ables computed with our new code and those computed 
with other codes is > 15-30%. In particular, the dis- 
agreement with the models computed by Young (2005) 
is probably due to the initial density profile. Indeed, in 
this investigation we limited our analysis to simplified ini- 
tial conditions, assuming that gas behaves as a polytrope 
and that p cx cx [sm(7ra;)/(7rx)]'^/* (see Eq. [6]). How- 
ever, the external layers of a SN typically have a steeper 
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Fig. 4. — Evoluti on of the main observables for model (1) in Table 
Eland model (A) of lYound ||2004I 1 (see text for details). Top, middle, 
and bottom panels show the bolometric light curve, photospheric 
velocity and temperature as a function of time, respectively. 




Fig. 7. — Same as Fig.|4l but for the comparison between model 
(10) and the corresponding model computed with the semi-analytic 
code (see text for details). 
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Fig. 5. — Same as Fig. [J] but for the comparison between model 
(3) and the corresponding model computed with the semi-analytic 
code (see text for details). 
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Fig. 8. — Same as Fig. [J] but for the comparison between model 
(13) and the corresponding model computed with the semi-analytic 
code (see text for details). 




0, 20000 

= 15000 

2 10000 

g. 5000 

E 



50 



100 150 
time [day] 




mod 4 

semi-analytic modelling 



100 150 
time [day] 



mod 4 

semi-analytic modelling 



100 150 
time [day] 



200 



250 



Fig. 6. — Same as Fig. |4l but for the comparison between model 
(4) and the corresponding model computed with the semi-analytic 
code (see text for details). 



power-law distribution caused by the acceleration of the 



1e-^43 t 

1e-^42 

1e-^41 

1e-^40 





5000 



20000 
15000 
10000 
5000 




mod 13 
Young's data 



50 



100 



150 200 
time [day] 



250 



300 



350 




40 60 

time [day] 



mod 13 
Young's data 



20 



40 



60 

time [day] 



80 



100 



120 



Fig. 9. — Same as Fig. |4] but for the co mparison between our 
model (13) and model (B) of lYounel l l200l) (see text for details). 



shock wave through the exponentially decaying density 
profile of the progenitor star. These layers (0.01-0.1 Mq 
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Fig. 10. — Bolomctric light curves of model s (1) and (2) compared 
to the lig ht curve of SN 1987A (taken from lCatchpole et al inMTl. 

mand IHamuv et aLllTOSST) . The best fitting model of lYounel 
) is also reported for comparison. 
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Fig. 11. — Evolution of the photospheric velocity (top panel) and 
temperature (bottom panel) for models (1) and (2) compared to the 
corresponding qu antities of SN 198 7A (taken from j Catchpole et al.l 
[1987, 1988, Ham uv et aLllTMa and lPhillips et allll988l) . 

in mass) are set into high -velocity and are not homolo - 
gously expanding (see also lUtrobinll2004l : IWooslevlll988D . 
During the first ~10-20 days after the explosion, the pho- 
tosphere is located in this shell, which is then cooler and 
has a higher velocity than the bulk of the ejecta. This 
results in a different evolution of the photospheric tem- 
perature and velocity. A better agreement during this 
phase is found with the semi-analytic code, that adopts 
a uniform density distribution. Some disagreement is 
present also at late-times (> 80-110 days), that is prob- 
ably related to differences in the treatment of ionization 
balance and in the adopted opacities (including the opac- 
ity floor). 

Despite the limitations of our simplified initial condi- 
tions, also the comparison with the bolomctric light curve 
of SN 1987A is satisfactory. We are able to reproduce its 
main features (peak luminosity and phase at maximum) 
with models having initial radius of 3 x 10^^ cm, total ini- 
tial energy of 1 foe, amount of ^^Ni equal to 0.07Mq, and 
envelope mass ranging between 16 and 18 Mq (models 
(1) and (2) in Figures [TOl and fTTj) . The agreement can 



be considered satisfactory also because we did not per- 
form any "fine-tuning" of the initial composition which 
is typically needed in order to accurately reproduce the 
observed shape (width and rise t o peak) of the bolomct- 
ric lig ht curve of SN 1987A (e.g. lWoosley"1988l; lUtrobh] 
l2Q04f l . Further residual differences may also be caused by 
the absence of non-thermal ionization from gamma rays 
in our models. The time evolution of the photospheric 
velocity and temperature of SN 1987A is also well repro- 
duced by models (1) and (2), apart from the differences 
in the early and the late-time evolution due to the same 
reasons mentioned above. Moreover, as can be seen from 
Figure [TUl the luminosity in the radioactive tail predicted 
by the model is lower than the observed one. This is a 
consequence of fallback occurring during the evolution. 
We found that the innermost ^ O.OIM0 of the envelope, 
containing ^ 2.4 x lO~^Af0 of ^^Ni, have been accreted 
onto the central remnant. 

4.2. Reference case 

To give a general overview of the post-explosion evo- 
lution of a "typical" CC-SN, we focus on the properties 
of our model (13) that has rather common initial pa- 
rameters. The evolution is determined by the thermody- 
namics of the expanding ejecta. The internal energy de- 
posited by the shock wave and that released by gamma- 
ray radioactive decays are used to expand the ejecta 
and power light curve. The evolution is characterized 
by three phases in which different heating and emission 
mechanisms dominate. During the first phase (diffusive 
phase), the envelope is completely ionized and optically 
thick, and the emission is due to the release of internal 
energy on a diffusion timescale. In the second phase (re- 
combination phase), the ejecta are recombining and the 
emission is dominated by the sudden release of energy 
caused by the receding motion of the wavefront through 
the envelope. During the last phase (radioactive-decay 
phase or radioactive tail), the envelope is recombined 
and optically thin to optical photons, and the emission 
comes from the thermalization of the energy deposited 
by gamma-ray photons. Observationally, the first two 
phases coincide with what is usually defined the plateau 
or photopheric phase, while the last phase is referred to 
as nebular phase. 

Figures [12] through [14] show the physical properties of 
model (13) at three different times, taken to be repre- 
sentative of the aforementioned three phases. Moreover 
Figures [15] through [18] show the evolution of the photo- 
pheric radius, and the radial profiles of the photospheric 
velocity and temperature in more detail. 

During the first phase (< 35 days), the radiation diffu- 
sion time-scale is much longer than the expansion time- 
scale, and the cooling induced by photon diffusion is neg- 
ligible. Although the internal energy decreases because 
of expansion, the temperature and the density are suffi- 
ciently high that the envelope remains completely ionized 
and optically thick (Figures [T2l [15] and [TB]) . The photo- 
pheric radius, which is located in the outermost shell, 
increases by a factor ^ 40 in ^ 35 days from the break- 
out of the shock wave at the stellar surface (see Figure 
lisp . The expansion is nearly homologous, as witnessed 
by the nearly constant value of the photospheric velocity. 

At 35 days, in the outermost layers the temperature 
reaches ~ 6000 K and hydrogen starts to recombine. This 
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Fig. 12. — From bottom to top, radial profiles of tlie gas matter 
density (in units of g cm~^), absolute value of the fluid radial- 
velocity (in units of c), radiative fiux (in units of erg cm~^), gas 
temperature (in units of Kelvin), and radiation energy density (in 
units of erg cm^'^) as a function of radius (in units of cm) . All scales 
on the y axis are logarithmic. The radial profiles refer to model 
(13) at 18 days from the breakout of the shock wave at the stellar 
surface. The innermost part of the envelope, below logii ~ 12.7, 
is accreting onto the central remnant and the radial-velocity there 
is negative. 
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Fig. 14. — Same as Fig. [12] but for 130 days from the breakout 
of the shock wave at the stellar surface. The boundary of the 
innermost accreting part of the envelope is at log R ~ 13.4. 




Fig. 15. — Gas temperature as a function of radius for model 
(13) at different times. Labels indicate the time (in units of days) 
from the breakout of the shock wave at the stellar surface. 
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Fig. 13. — Same as Fig. 1121 but for 75 days from the breakout 
of the shock wave at the stellar surface. The boundary of the 
innermost accreting part of the envelope is at log/? ~ 13.3. 

marks the beginning of the recombination phase. Dur- 
ing this phase, the evolution is characterized by the fast 
recombination of the outer zone, leading to the forma- 
tion of a recombination wave (RW) that moves inward 
(in mass) like a flame. In the external layers, the tem- 
perature drops below 6000 K and, consequently, the 
internal energy is efficiently radiated away because of the 
sudden decrease of the opacity (see Figures [131 [T5] and 
[TB)) . The RW marks the boundary between the inner en- 
velope that is optically thick, ionized and hot (> 6000 K), 
and the outer layers that are optically thin, recombined 
and cooler (< 6000 K). The photosphere moves inward 
(in mass) following the motion of the RW and allows pho- 
tons to escape sooner than they would if the photosphere 
were fixed in the outmost layer of the envelope. As a con- 
sequence, advective cooling induced by the WR motion 
becomes dominant. The energy radiated away at this 





5d 






lOd 






15d -- 






18d 






. 25 d 






" -.^ 35 d 




_ 


bOd 






- 75 d 






i 'r " ^- ■'llOd 






:130d 











5 10 15 20 

M [M3„„l 

Fig. 16. — Same as Fig. 1151 but for the gas temperature as a 
function of mass. 

stage is the sum of the residual internal energy left over 
after expansion and that liberated during recombination. 
The energy deposited by gamma-rays becomes signifi- 
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Fig. 17. — Absolute value of the fluid radial-velocity as a function 
of radius for model (13) at different times. Labels indicate the time 
(in unit of days) from shock breakout. The minimum of each curve 
marks the boundary of the inner accreting zone (where the radial- 
velocity is negative). 
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Fig. 18. — Evolution of the photospheric radius for model (13). 



cant only when t he RW reaches the innermost zon es, full 
of ^^^Ni (see also lYounel[20p : iBersten et al.|[20TTI ). The 
bolometric luminosity, initially almost constant, plum- 
mets when the RW reaches the center (see the top panel 
in Figure in]). The bulk of the envelope (^ 95% in mass) 
is recombined after ~ 110 days from the breakout of the 
shock wave at the stellar surface (see Figure [T5|) . The 
photopheric velocity decreases because the photosphere 
moves inward with the RW, reaching inner and slower 
layers as time elapses (see the middle panel in Figure IH]). 

The velocity of propagation of the RW is determined 
by the physical conditions in the envelope, in particular 
the temperature, density, amount of ^^Ni and their dis- 
tributions. The motion of the photospheric radius for a 
distant observer is the result of the interplay between the 
inward motion of the RW in mass and the expansion of 
the ejecta. Consequently, as shown in Figure [T51 from ^ 
35 up to ~ 65 days, the photospheric radius increases, 
while later on it stays at an almost constant (eulerian) 
radial coordinate for about ~ 20 days. Then, the pho- 
tosphere starts to recede in radius and the luminosity 
plummets. 



The duration and steepness of the transition from the 
plateau t o the radioactive tail depend on the ^^Ni mass 
(see sect. 14. 3p and its distribution in the ejecta. During 
this last phase, the evolution is completely governed by 
the deposition and reprocessing of the energy released 
from the radioactive decay of ^^Co into ^^Fe. The en- 
velope is sufficiently opaque to the gamma-rays emitted 
in the radioactive decays that a large fraction of them 
are effectively absorbed and their energy deposited in 
the ejecta. This energy is easily radiated away by lower 
energy photons on a timescale much shorter than the 
expansion timescale. Also the heating time, which is 
essentially the ^^Co decay time, is smaller than the ex- 
pansion time and, hence, the loss of internal energy due 
to expansion is negligible. Therefore, the bolometric lu- 
minosity is the total heating rate from the radioactive 
decays, proportional to the mass of ^^Ni synthesized in 
the explosion (see also sect. 



4.3. Effects of varying the initial parameters on the 
light curve, photospheric velocity and temperature 

The sample of models (reported in Table ^ enabled 
us to study the effects of the variation of each param- 
eter (namely, envelope mass Menv, initial outer radius 
i?o, total ejecta energy E, and amount of ^^Ni Mmi) on 
the main observables, while holding the others fixed. Fig- 
ures[Tn]through[21]summarize the results of this analysis. 
We note, however, that this investigation is not meant 
to be a realistic model survey for a detailed quantita- 
tive comparison with observations, because it does not 
include models computed starting from realistic initial 
conditions. This is postponed to a follow-up investiga- 
tion. 
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Fig. 19. — Effects of varying the envelope mass on the main 
observables for models (3), (4), and (10) with Menv = 8, 12, and 
16Mq , respectively. The models (see Tabled have the same initial 
radius (i?o = 3 X 10^^ cm), total energy (E = 1 foe), and amount 
of ^'^Ni (Mivi = O.O35M0). Top, middle, and bottom panels show 
the bolometric light curve, photospheric velocity and temperature 
as a function of time, respectively. 



As can be seen from Figure [TOl more massive envelopes 
lead to a fainter bolometric light curve during the diffu- 
sive phase, a longer plateau phase, a hotter photospheric 
temperature during the first ~ 20-30 days, a lower pho- 
tospheric velocity up to 50 days and a slower decline at 
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Fig. 20. — Same as Fig. \W\ but for the initial radius. The lines 
refer to models (1), (13), (21), and (22) with i?,o = SxlO^^, 3x10^^, 
6 X 10^'^, and 10^* cm, respectively. The models (see Table[2)l have 
the same envelope mass (Menv = IGMq), total energy {E = I 
foe), and amount of S^Ni (Mjvi = O.OTM©). 
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Fig. 21. — As for Fig. ll9l but for the total energy. The lines refer 
to models (6), (13), (14), (15), (16), (18), and (20) with E = 0.5, 
1, 1.3, 1.6, 1.9, 2, and 3 foe, respectively. The models (see Table 
[2)1 have the same initial radius {Rq = 3 x 10^^ cm), envelope mass 
{Menv = I6M0), and amount of s^Ni (Mffi = O.OTMq). 

longer times. As known, this behavior is a consequence 
of the increase in the diffusion timescale with increasing 
envelope mass, that causes a slower leakage of photons 
and hence a lower luminosity, longer plateau and higher 
temperature. At the same time, increasing the mass and 
keeping the initial kinetic energy constant (which is half 
of the total energy; see Sect. I3.2p results in a smaller 
photospheric velocity during the first weeks and a slower 
decline later on, caused by the slower motion of the RW 
in a higher density envelope. 

The effects of varying the initial radius are shown in 
Figure [20I Larger values of Rq cause a higher luminosity 
during the diffusive and recombination phases, a longer 
plateau, and a smaller photospheric velocity up to ~ 25 
days. The higher luminosity L during the diffusive phase 
is caused by the fa ct that L is p roportional to the ini- 
tial radius Rq (e.g. lArnettI [19960 . as the radiative flux 
F (X T'^ /{kssPoRo) oc Rq . The lower initial photo- 
spheric velocity, longer plateau and higher initial temper- 
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Fig. 22.— As for Fig.[l9l but for the ^^Ni mass. The lines refer to 
models (7), (8), (9), (10), (11), (12), and (13) with Mjvi = 0.005, 
0.01, 0.02, 0.035, 0.045, 0.055, and O.OTOM©, respectively The 
models (see Tabled have the same initial radius (Ho = 3 X 10^^ 
cm), envelope mass (Me„„ = I6M0), and total energy {E = I 
foe). 

ature are a consequence of the smaller PdV work done 
by bigger expanding envelopes. Less internal energy is 
initially converted into kinetic energy, the expansion is 
slower and more energy remains available for heating up 
the envelope. 

Figure shows the effects of varying the total energy 
(with ratio between thermal and kinetic energy equal to 
1; see Sect. 13.21) . Increasing it causes a higher luminos- 
ity in the bolometric light curve during the diffusive and 
recombination phases, a shorter plateau, a higher photo- 
spheric velocity up to ~ 60 — 70 days and a faster decline 
at longer times. The early behavior can be easily un- 
derstood as a consequence of the larger internal energy 
dumped in the ejecta, that increases their initial velocity 
and accelerates all the evolutionary stages. 

As for the effects of varying ^^Ni (see Figure . in- 
creasing it makes the plateau phase longer and causes 
a slower decline of the bolometric light curve towards 
the radioactive tail. Furthermore, as Mjvi increases, also 
the decline of the photospheric velocity is slower, while 
the bolometric luminosity becomes higher during the ra- 
dioactive tail. All these effects are related to the in- 
creased heating provided by the larger amount of ^^Ni 
concentrated in the innermost part of the envelope, that 
increases the internal energy and luminosity especially at 
late phases and slows down the motion of th e RW (e.g. 
lArnettI [19961 IHamuvl [20031 INadvozhh3[2003l) . 

5. SUMMARY AND FINAL COMMENTS 

We developed a general-relativistic, radiation hydro- 
dynamics, Lagrangian code tailored to the radiation- 
hydrodynamical modelling of CC-SNe, whose distinctive 
features are an accurate treatment of radiative transfer 
coupled to relativistic hydrodynamics, a self-consistent 
treatment of the evolution of the innermost ejecta tak- 
ing into account the gravitational effects of the central 
compact remnant, and a fully implicit Lagrangian ap- 
proach to the solution of the coupled non-linear finite 
difference system of equations. The latter property rep- 
resents the major novelty of this code, that is based on 
a previous version originally developed for studying fall 
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back in the aftermath of a SN explosion ([Zampieri et alj 
[l998t IBalberg et alJ[2000l ). 

With this code, a total of 22 models of hydrogen-rich 
CC-SNe were calculated, enabling us to (a) validate the 
code against observations and similar models reported 
in the literature, and (b) study the role of the "main" 
parameters affecting the post-explosion evolution of the 
CC-SN events (namely the ejected mass, the progenitor 
radius, the explosion energy, and the amount of ^^Ni). 

The aforementioned grid of models permitted us also 
to explore possible correlations among different quan- 
tities that can be measured from the light curve and 
spectra, and investigate how variations of basic pa- 
rameters combine together to produce them. This 
is particularly relevant in connection with the pos- 
sibility to calibrate hydrogen-rich CC-SNe to turn 
them into usable distance indicators. Indeed, in re- 
cent years growing attention has been devoted to the 
construction of Hu bble diagrams using hydrogen- rich 
CC-S N events (e.g. iNugent et all 120061 : lOlivares et all 
120101 ) in order to derive cosmological parameters in- 
dependently of the usual method based on ty pe la 
SNe fe.g . i Riess et al.l 119981: iPerlmutter et all Il999l: 



Asticr et al 



20061: 1 W ood- V asev et alJl2007HHicken et al.l 



200a: .Freedman et al..i2009i) . 



Although a complete study using a more extended 
grid of models calculated from realistic initial conditions 
is presently under investigation, here we mention some 
preliminary results obtained for the sample of 22 mod- 
els presented in this paper. Indeed, despite the limita- 
tions related to the use of simplified initial conditions, 
we are able to reproduce the power-law relation between 
the luminosity and the photospheric velocity (both mea- 
sured at day 50 from the breakout of the shock wave 
at the stellar surface) foun d in the observational sample 
of iHamuv fc Pintol (|2002[ ). The index of the theoreti- 
cal correlation is equal to 3.13 ± 0.29, in good agree- 
ment withi n the errors wit h the value 3.03 ± 0.37 re- 
ported by IHamuv fc Pintol (|2002[ ). Our models con- 
firm also the anti-correlation between the light curve 
slope at the so-called inflection time ti (when the semi- 



logarithmic derivative of the luminosity L at the end of 
the plate au is maximum) and the amounts of ^ ^Ni in- 
ferred by lElmhamdi. Chugai. fc Danzigeii ()2003[ ) on ob- 
servational bases. Finally, we found a very promising cal- 
ibration relation between the luminosity L» at a generic 
time during the plateau and the characteristic time 
tc = io.4 — t*, where to. 4 is the time when decreases 
by a factor 2.5. These last two relations may represent 
useful tools for calibrating hydrogen-rich CC-SNe using 
only photometric data. In particular, for type II plateau 
SNe the latter relation is essentially independent of the 
explosion epoch as is approximately constant during 
the plateau. 

We are working at present to check the validity of 
these preliminary results against a more extended grid 
of models that is being computed from realistic initial 
conditions. We plan to check these calibrations also 
against observations using well sampled light curves of 
hydrogen-rich CC-SNe that are being collected within 
the dedicated large program "Supernova Variety and Nu- 
cleosynthesis Yields" (P.I. S. Benetti) presently running 
at the European Southern Observatory and the Telesco- 
pio Nazionale Galileo. The aforementioned grid of mod- 
els will also serve to build an extended database to be 
compared with observations of single SNe in order to in- 
fer their physical properties, in analogy to what already 
being done using models with simphfied initial conditions 
(e.g. SN 2007od; see In serra et al1l2011| ). 

Our medium- and long-term goal is the development of 
a sort of "CC-SNe Laboratory" in which our code is inter- 
faced, in input, with other codes deali ng with the calcula- 
tions of the pre-SN evolution (see e.g. iLimongi fc Chieffil 
120031 and refer ences therein) and explo sive nucleosyn- 
thesis (see e.g. iChieffi fc Limongil 12004 and references 
therein), and in output with a spectral sy nthesis code 
(see e.g. lMazzalill2000l : iMazzah fc L"iIcvl [T993). This will 
allow us to describe the evolution of a CC-SN event in 
a "self-consistent" way from the evolutionary stages pre- 
ceding the main sequence up to the post-explosive phases 
as a function of initial mass, metallicity, stellar rotation, 
and mass loss history of the CC-SN progenitor. 



APPENDIX 
FINITE DIFFERENCE EQUATIONS 

We start from Eqs. [H , |2| and |3l , write r_t = au, and express = —{pr'^)^t/pT'^ in terms of the continuity equation 
(Eqs. [8] and [4] iu iZampieri et al1ll998D . We then discretize the spatial computational domain (in Lagrangian mass 
ji) dividing it into a grid of jmax zones (equal to 110 in our simulations), where a constant fractional increment in grid 
spacing between successive zones a is used. This fractional increment is calculated from the following relation: 



Mr) 



E 



(Al) 



where fimin — P-jmin is the inner boundary of the grid in Lagrangian mass, Umax = Pj^a^ is the outer boundary (fixed 
during a simulation due to the conservation of the envelope mass), and Ap,j_^_i/2 = fJ-j+i — fJ'j — At the 

beginning of the simulation, fimin = 0, Pmax — ^env, and the mass contained within the first shell A/ii/2 is fixed by 
the requirement that the radial spacing between the first two shells is 30% of the inner radial boundary ri„ = r{pmin)- 
During the evolution, if the inner edge of the innermost zone crosses ri„, it is removed from the calculation and a 
regridding of all the variables is performed so as to have jmax zones at all times. In particular, the new inner boundary 
is set equal to pj^^^-^-i^ where / is the number of zones that have crossed r^, and the new fractional increment a is 
calculated from Eq. |Alj with the new value of inner boundary. Afterwards, all the variables are interpolated on the 
new grid. 

As for the spatial centering of the variables, B — auT^ and wq are evaluated at mid zone {pj-1/2), while wi at 
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zone boundary (/ij). Concerning the time centering, all quantities are evaluated at the full time level t". With this 
centering and following a standard Lagrangian approach for the discretization of the equations and the derivatives, 
the finite difference form of Eqs. [T], [2] and [3] can be written as: 



n+l _ n 

Ai"+i/2 



(^0);-l/2 - (^0)"-l/2 



' a 



-1/2 
J-l/2 



+ 1/2 
J-1/2 



n+l/2 
^3-1/2 
At"+V2 



n+l/2 



ri+l 
J-1/2 



j-1/2 ('«0)J_i/2 



0-1/2, 



-1/2 ^'^-P''j-l/2Pj-l/2 



B 



n+l/2 
i-1/2 



1+1/2, 



-1/2 



j-l/2"'j-l/2 



+ l/2^^n+l/2^2 



1/2 

J-1/2 



1 



n+l/2 
1/2 



,'J,2^"-^l/2 



"-H/2/ "-M/2N2 
7-1 I 



-1/2 



-1/2 



'J~l 



_ n+l/2 

r/j-1/2 



2+1/2 r 



J-1/2 
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-1/2 



0-1 



Airra ( 4 



-1 n+l/2 



= 



J-l/2 



^^n+1/2 



^^n+1/2 



Ai"+l/2 



= 0, 



(A2) 



(A3) 
(A4) 



where 



B 



Tya I r>n+l 

1+1/2 _ -"^-1/2 +^j-l/2 

J-1/2 2 



(^^)n+l/2_(-0)"-l/2 + K);+V2 



''j-1/2 2 

F,„=Ar+i/2a"+i/2 



2 Mj + l/2?'j + l/2 - ■"j-l/2''i-l/2 



n+l/2 



B,„ = At 



n+l/2 1 87raj-rj- 



7'j + l/2 - ''i-l/2 
(1/3 + /j + 1/2) (wo)j + l/2 - (1/3 -f /j-1/2) (wo)j-l/2 



/ '^i+1/2 - ''j-1/2 



rj+1/2 - rj-1/2 

Hl/2 



(A5) 
(A6) 
(A7) 



(A8) 



As £, jp, A:j^ and kp are in general rather complex, non- linear functions of p and B (see, e.g., Eqs. [14] and [15] 
in lZampieri et al.lll998[ ). relations |A2| - |A4| form a highly non-linear system of (Sjmax- ^) equations havi ng (3jmaa;-3) 
unknowns, that we solve by means of the Newton-Raphson iterative method (see, e.g.. iPress et "aI1 [l996l for details). 
Using this approach, the finite difference equations are linearized and solved iteratively for the corrections at the new 
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time level starting from an initial guess. The equations for the corrections are: 
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Ai3,-i/2 + -3^^^f^AK),-l/2 



Ai3,-i/2 + -g|e^^AK),-i/2 



c'('i«i ). 

^«F^AK). 



^Jniin + 1) 



[2](J) 
[3](J) 



where J j'^^ j), i j'^^ j), and i j'^^ j) are the left-hand sides of the Eqs. |A2| . |A3| . and |A4| . respectively and the value 
of their derivatives is calculated numerically. This system of equations can be written in matrix form as: 



Ax 



(A9) 



where A is the matrix of the coefhcients (the so-called Jacobian of the system), the raised dot denotes matrix mul- 
tiplication, and X and b are the unknown corrections and the known right-hand sides (written as column vectors), 
respectively. They read: 
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AB( 



-1/2 



^(^0)0,„.„ + l)-l/2 

A(w;o)j_i/2 



A(wi)(j- 



-1/2 

.ox) 
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Given the functional dependence of i j'^^ j), i j'^^ j), and i j'^^ j) from i?j-i/2, (u'o)j"-i/2i and (wi)j, the matrix 
of coefficients A is a so-caUed "band diagonal" matrix having nonzero elements only on the diagonal plus mi (=6 at 
most) elements immediately to the left of (or below) the diagonal and m2(=4 at most) elements immediately to the 
right (or above it). In other words, the elements of A are equal to zero if either the column index h is greater than 
the row index k + TO2 or the row index k is greater than the column index h + mi . 

In order to solve the system |A9| in an efficient way, minimizing the CPU time and the required storage space, we 
adopt a matrix inversion package based on a LU decomposition, where the matrix A is stored an d manipulated in a 
so-called "compact form", in which only the non-zero elements are considered (see, e.g.. Press et al. 1996!', for details). 

In addition, in order to increase the numerical accuracy, we use an iterative improvement of the solution x. Specif- 
ically, we approximate x with Xq -I- i5xo, where Xq is the solution of the syst em A ■ x = b and (5xo - the first-order 
correction to Xq - is the solution of the system A • (5x = (b — Axq) (see, e.g.. lPress et al.|[l996l for details). 
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